home *** CD-ROM | disk | FTP | other *** search
/ PD Collection CD 1 / PD Collection CD 1.iso / programer2 / pari2 / pari / c / polarit2 < prev    next >
Text File  |  1991-11-28  |  41KB  |  1,439 lines
  1. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*~                                                                     ~*/
  4. /*~                       OPERATIONS ARITHMETIQUES                      ~*/
  5. /*~                                                                     ~*/
  6. /*~                           SUR LES POLYNOMES                         ~*/
  7. /*~                                                                     ~*/
  8. /*~                           (deuxieme partie)                         ~*/
  9. /*~                                                                     ~*/
  10. /*~                          copyright Babe Cool                        ~*/
  11. /*~                                                                     ~*/
  12. /*~                                                                     ~*/
  13. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  14. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  15.  
  16. # include "genpari.h"
  17.  
  18. /* setup for calling cmbf */
  19.  
  20. GEN cmbf_target, cmbf_lc, cmbf_abslc, cmbf_abslcxtarget, cmbf_mod,
  21.     cmbf_modfax, cmbf_fax;
  22. long cmbf_degree, cmbf_modfaxn, cmbf_faxn;
  23.  
  24. #define NOMBDEP 5
  25.  
  26. GEN releve(x)
  27.      GEN x;
  28. {
  29.   long lx=lg(x),tx=typ(x),i;
  30.   GEN y;
  31.   switch(tx)
  32.     {
  33.     case 3:
  34.     case 9: return (GEN)x[2];
  35.     case 7: return (GEN)x[4];
  36.     case 10: lx=lgef(x);
  37.     case 11: if(!signe(x)) return x;
  38.     case 6:
  39.     case 17:
  40.     case 18:
  41.     case 19: y=cgetg(lx,tx);
  42.       for(i=1;i<lontyp[tx];i++) y[i]=x[i];
  43.       for(i=lontyp[tx];i<lx;i++)
  44.     y[i]=(long)releve(x[i]);
  45.       return y;
  46.     default: return x;
  47.     }
  48. }
  49.  
  50. GEN centermod(x,p)
  51.      GEN x,p;
  52.  
  53. {
  54.   GEN y,p1,ps2;
  55.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i;
  56.  
  57.   ps2=shifti(p,-1);
  58.   switch(tx)
  59.     {
  60.     case 1: y=modii(x,p);
  61.       if(gcmp(y,ps2)>=0)
  62.     {tetpil=avma;y=gerepile(av,tetpil,subii(y,p));}
  63.       return y;
  64.     case 10: lx=lgef(x);
  65.     case 11: if(!signe(x)) return x;
  66.       y=cgetg(lx,tx);
  67.       y[1]=x[1];for(i=2;i<lx;i++)
  68.     {
  69.       p1=modii(x[i],p);
  70.       if(gcmp(p1,ps2)>=0) p1=subii(p1,p);
  71.       y[i]=(long)p1;
  72.     }
  73.       tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  74.     case 6:
  75.     case 17:
  76.     case 18:
  77.     case 19: y=cgetg(lx,tx);
  78.       for(i=1;i<lx;i++) y[i]=(long)centermod(x[i],p);
  79.       return y;
  80.     default: return x;
  81.     }
  82. }
  83.  
  84.  
  85. /* This code was kindly written for us by Richard Schroeppel */
  86.  
  87.  
  88. /* Note that PARI's idea of the maximum possible coefficient involves the
  89. limit on the degree (klim).  Consider revising this.
  90. If I don't respect the degree limit when testing potential factors,
  91. there's the possibility that I might identify a high degree factor that
  92. isn't irreducible, because it's lower degree divisors were missed because
  93. they had a coefficient outside the Borne limit for klim, but the higher
  94. degree factor had it's coefficients within Borne.  This would still have
  95. the property that any factors of degree <= klim were guaranteed irr, but
  96. higher degrees (> 2*klim) might not be irr. */
  97.  
  98.  
  99.  
  100. /* the subroutine:
  101.    fxn points at the first unconsidered factor for the current combination
  102.    psf is the product-so-far, or 0 for a null product
  103.    dlim is the degree limit remaining for unconsidered divisors
  104.    other arguments are "global" and must already be setup
  105.    as factors are found, they are put in cmbf_fax; the count is kept in
  106.    cmbf_faxn; and they are divided out of cmbf_target; the degree and
  107.    leading coefficient are updated; and the constituent modular factors
  108.    are deleted from cmbf_modfax.
  109.    exit value is 1 if any factors are found.
  110.    If psf is 0, all factors made up from pieces at or after fxn will be
  111.    found & removed.  If psf is not 0, only the factor which is a
  112.    continuation of psf will be found.
  113.   */
  114.  
  115. int combine_factors(fxn,psf,dlim)
  116.      long fxn,dlim;
  117.      GEN psf;
  118.  
  119. {
  120.   int val=0, val2=0;  /* assume failure */
  121.   GEN newf, newpsf, quo, rem, cont; long newd;
  122.  
  123.   if (dlim <= 0) return 0;
  124.   if (fxn > cmbf_modfaxn) return 0;
  125.   /* first, try deeper factors without considering the current one */
  126.   if (fxn < cmbf_modfaxn) { val = combine_factors(fxn+1,psf,dlim);
  127.                 if (val&&psf) return val; };
  128.   /* second, try including the current modular factor in the product */
  129.   newf = (GEN)cmbf_modfax[fxn];
  130.   if (!newf) return val;  /* modular factor already used */
  131.   newd = lgef(newf)-3;
  132.   if (newd > dlim) return val;  /* degree of new factor is too large */
  133.   if (psf) { newpsf = centermod(gmul(psf,newf),cmbf_mod); }
  134.      else { newpsf = centermod(gmul(cmbf_abslc,newf),cmbf_mod); };
  135.   /* try out the new combination */
  136.   quo = poldivres(cmbf_abslcxtarget,newpsf,&rem);
  137.   if (!signe(rem))  /* found a factor */
  138.    { 
  139.      cont = content(newpsf);  /* hope this is always positive */
  140.      newpsf = gdiv(newpsf,cont);
  141.      cmbf_fax[++cmbf_faxn] = (long)newpsf; /* store factor */
  142.      cmbf_modfax[fxn] = 0; /* remove used modular factor */
  143.      /* fix up target */
  144.      cmbf_target = gdiv(quo,newpsf[lgef(newpsf)-1]);
  145.      cmbf_degree = lgef(cmbf_target)-3;
  146.      cmbf_lc = (GEN)cmbf_target[cmbf_degree+2];  /* leading coefficient */
  147.      cmbf_abslc = gabs(cmbf_lc); /* |lc| */
  148.      cmbf_abslcxtarget = gmul(cmbf_abslc,cmbf_target); /* abslc * target */
  149.      return 1;
  150.    }
  151.   /* newpsf needs more; try for it */
  152.   if (newd == dlim) return val;  /* no more room in degree limit */
  153.   if (fxn == cmbf_modfaxn) return val;  /* no more modular factors to try */
  154.   val2 = combine_factors(fxn+1,newpsf,dlim-newd);
  155.   if (val2) cmbf_modfax[fxn] = 0; /* remove used modular factor */
  156.   return val||val2;
  157. }
  158.  
  159.  
  160. GEN  squff(a,klim)
  161.      GEN a;
  162.      long klim;
  163. {
  164.   GEN borne,di,p1,p2,y,g,pe,pg,s,t,u,v,w,xmod,unmod,polxmodp;
  165.   GEN ae,be,aprov,pev;
  166.   GEN tabd[NOMBDEP][500],unmodp[NOMBDEP];
  167.   unsigned long tabbit[60],tabkbit[60],tablbit[60],j1,rem,pro;
  168.   byteptr pt=diffptr;
  169.   long av=avma,tetpil,p=0,tabp[NOMBDEP],nfacp[NOMBDEP];
  170.   long va=varn(a),da=lgef(a)-3,lbit,k,d0,fla,i,j,d,e,d1,d2;
  171.   long np,nmax,imax,kd,lgg,nf,ev,nfd,nft;
  172.   
  173.   di=discsr(a);np=0;lbit=(da>>4)+1;nmax=da+1;imax=0;kd=da>>1;
  174.   if((!klim)||(klim>kd)) klim=kd;
  175.   while(np<NOMBDEP)
  176.     {
  177.       p+= *pt++;
  178.       if(signe(gmodgs(di,p))&&signe(gmodgs(a[da+2],p)))
  179.      /* require p doesn't divide leading coefficient */
  180.     {
  181.       tabp[np]=p;unmodp[np++]=gmodulcp(un,stoi(p));
  182.     }
  183.     }
  184.   for(i=0;i<NOMBDEP;i++)
  185.     {
  186.       p=tabp[i];nfacp[i]=0;unmod=unmodp[i];
  187.       tabkbit[lbit-1]=1;for(j=0;j<lbit-1;j++) tabkbit[j]=0;
  188.       d=0;v=gmul(unmod,a);xmod=gmodulcp(polxmodp=gmul(unmod,polx[va]),v);w=xmod;
  189.       while(d<(e=(lgef(v)-3))>>1)
  190.     {
  191.       d++;
  192.       w=gpuigs(w,p);p1=(GEN)(gsub(w,xmod))[2];
  193.       g=ggcd(p1,v);
  194.       tabd[i][d]=g;lgg=lgef(g)-3;
  195.       if(lgg>0)
  196.         {
  197.           v=gdiv(v,g);w=gmodulcp(gmod(w[2],v),v);
  198.           xmod=gmodulcp(polxmodp,v);
  199.           nfacp[i]+=(lgg/d);
  200.           for(kd=d;kd<=lgg;kd+=d)
  201.         {
  202.           d1=d>>4;d2=d&15;rem=0;
  203.           for(j=0;j<d1;j++) tablbit[lbit-1-j]=0;
  204.           for(j=d1;j<lbit;j++)
  205.             {
  206.               pro=tabkbit[lbit-1-j+d1]<<d2;
  207.               tablbit[lbit-1-j]=(pro&0xffff)+rem;rem=pro>>16;
  208.             }
  209.           for(j=0;j<lbit;j++) tabkbit[j] |= tablbit[j];
  210.         }
  211.         }
  212.     }
  213.       if(e>0)
  214.     {
  215.       tabd[i][e]=v;nfacp[i]++;
  216.       d1=e>>4;d2=e&15;rem=0;
  217.       for(j=0;j<d1;j++) tablbit[lbit-1-j]=0;
  218.       for(j=d1;j<lbit;j++)
  219.         {
  220.           pro=tabkbit[lbit-1-j+d1]<<d2;
  221.           tablbit[lbit-1-j]=(pro&0xffff)+rem;rem=pro>>16;
  222.         }
  223.       for(j=0;j<lbit;j++) tabkbit[j] |= tablbit[j];
  224.     }
  225.       if(nfacp[i]<nmax) {nmax=nfacp[i];imax=i;}
  226.       for(j=d+1;j<e;j++) tabd[i][j]=polun[va];
  227.       if(i)
  228.     for(j=0;j<lbit;j++) tabbit[j] &= tabkbit[j];
  229.       else
  230.     {
  231.       for(j=0;j<lbit;j++) tabbit[j] = tabkbit[j];
  232.     }
  233.       fla=0;j1=1;
  234.       for(j=lbit-1;(j>=0)&&(!fla);j--)
  235.     {
  236.       for(k=0;(k<=15)&&(!fla);k++)
  237.         {
  238.           fla=tabbit[j]&j1;
  239.           if((!k)&&(j==lbit-1)) fla=0;
  240.           j1<<=1;
  241.         }
  242.       j1=1;
  243.     }
  244.       d0=((lbit-j-2)<<4)+k-1;
  245.       if((d0==da)||(d0>klim))
  246.     {
  247.       tetpil=avma;y=cgetg(2,18);y[1]=lcopy(a);
  248.       return gerepile(av,tetpil,y);
  249.     }
  250.     }
  251.   p1=gzero;
  252.   for(i=2;i<=da+2;i++) p1=gadd(p1,gmul(a[i],a[i]));
  253.   p1=gadd(p2=gabs(a[da+2],4),gaddsg(1,racine(p1)));
  254.   borne=gmul(p1,binome(stoi(da-1),klim));
  255.   p=tabp[imax];unmod=unmodp[imax];
  256.   e=itos(gceil(gdiv(glog(gmul2n(gmul(p2,borne),1),4),glog(pg=stoi(p),4))));
  257.   pev=gpuigs(pg,e);
  258.   y=cgetg((nf=nfacp[imax])+1,18);k=0;
  259.   p1=cgetg(nf+1,18);nft=0;
  260.   for(d=1;nft<nf;d++)
  261.     {
  262.       g=tabd[imax][d];g=gdiv(g,g[lgef(g)-1]);lgg=lgef(g)-3;
  263.       if(lgg)
  264.     {
  265.       nfd=lgg/d;p1[nft+1]=(long)g;
  266.       split(p,p1+nft+1,d,p,shifti(gpuigs(pg,d),-1));
  267.       nft+=nfd;
  268.     }
  269.     }
  270. /* do a Hensel lift */
  271.   aprov=a;
  272.   for(i=1;i<=nf-1;i++)
  273.     {
  274.       pe=pg;ev=1;p2=(GEN)p1[i];
  275.       ae=gdiv(p2,p2[lgef(p2)-1]);be=gdiv(aprov,ae);
  276.       g=bezoutpol(ae,be,&u,&v);
  277.       if(isnonscalar(g)) err(henser1);
  278.       u=gdiv(u,g[2]);v=gdiv(v,g[2]);
  279.       for(j=2;j<lgef(ae);j++) ae[j]=(long)releve(ae[j]);
  280.       for(j=2;j<lgef(be);j++) be[j]=(long)releve(be[j]);
  281.       for(j=2;j<lgef(u);j++) u[j]=(long)releve(u[j]);
  282.       for(j=2;j<lgef(v);j++) v[j]=(long)releve(v[j]);
  283.       do
  284.     {
  285.       g=gdiv(gsub(aprov,gmul(ae,be)),pe);
  286.       t=poldivres(gmul(v,g),ae,&s);for(j=2;j<lgef(s);j++) s[j]=lmod(s[j],pe);
  287.       ae=gadd(ae,gmul(pe,s));
  288.       s=gadd(gmul(u,g),gmul(t,be));for(j=2;j<lgef(s);j++) s[j]=lmod(s[j],pe);
  289.       be=gadd(be,gmul(pe,s));
  290.       g=gdiv(gsub(gun,gadd(gmul(u,ae),gmul(v,be))),pe);
  291.       t=poldivres(gmul(v,g),ae,&s);
  292.       for(j=2;j<lgef(s);j++)
  293.         s[j]=lmod(s[j],pe);
  294.       v=gadd(v,gmul(pe,s));
  295.       s=gadd(gmul(u,g),gmul(t,be));for(j=2;j<lgef(s);j++) s[j]=lmod(s[j],pe);
  296.       u=gadd(u,gmul(pe,s));
  297.       pe=gsqr(pe);ev<<=1;
  298.     }
  299.       while(ev<e);
  300.       p1[i]=(long)ae;if(i<nf-1) aprov=gdeuc(aprov,ae);
  301.     }
  302.   unmod=gmodulcp(gun,pev);
  303.   p1[nf]=(long)centermod(gmul(be,ginv(gmul(be[lgef(be)-1],unmod))[2]),pev);
  304.   for(i=1;i<=nf;i++) 
  305.     {g=(GEN)p1[i];if(!gcmp1(g[lgef(g)-1])) err(factpoler2);}
  306.   cmbf_target = a;  /* target poly.  Assumes content removed  */
  307.   cmbf_degree = lgef(cmbf_target)-3;
  308.   cmbf_lc = (GEN)cmbf_target[cmbf_degree+2];  /* leading coefficient */
  309.   cmbf_abslc = gabs(cmbf_lc); /* |lc| */
  310.   cmbf_abslcxtarget = gmul(cmbf_abslc,cmbf_target); /* abslc * target */
  311.   cmbf_mod = pev;   /* Modulus */
  312.   cmbf_modfax = p1;  /* array of modular factors.  Each has LC 1.
  313.             1 based indexing.  Product should be congruent to a/lc(a).  */
  314.   cmbf_modfaxn = nf;  /* number of modular factors */
  315.   cmbf_fax = cgetg(nf+2,18);
  316.   /* Result array.  Extra cell for leftover constant. */
  317.   cmbf_faxn = 0; /* pointer into result array; last used cell;
  318.             # of factors found */
  319.   
  320.   /* sorting factors decreasing by degree helps if klim is used */
  321.   /* if klim isn't used, can start with first arg of 2 instead of 1,
  322.      saving some time */
  323.   combine_factors(1,0,klim);  /* the call */
  324.  
  325.   /* follow-up */
  326.  
  327.   if (signe(cmbf_lc)<0) cmbf_target = gneg(cmbf_target);
  328.   if (cmbf_degree) cmbf_fax[++cmbf_faxn] = (long)cmbf_target; /* leftover factor */
  329. /*  if (signe(cmbf_lc)<0) cmbf_fax[++cmbf_faxn] = stoi(-1); */
  330.   tetpil=avma;y=cgetg(cmbf_faxn+1,18);
  331.   for(i=1;i<=cmbf_faxn;i++) y[i]=lcopy(cmbf_fax[i]);
  332.   return gerepile(av,tetpil,y);
  333. }
  334.   
  335. GEN factpol(x,klim)
  336.   GEN x;
  337.   long klim;
  338.  
  339. /* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs de degre <= klim */
  340.  
  341. {
  342.   long av=avma,av2,lx,vv,k,i,j,i1,f,nbfac;
  343.   GEN res,fa,p1,p2,y,d,a,ap,t,v,w;
  344.   
  345.   if((typ(x)!=10)||(!signe(x))) err(factpoler1);
  346.   y=cgetg(3,19);if((lx=lgef(x))==3) 
  347.     {y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  348.   if(lx==4)
  349.     {
  350.       p1=cgetg(2,18);y[1]=(long)p1;p1[1]=lcopy(x);
  351.       p1=cgetg(2,18);y[2]=(long)p1;p1[1]=un;return y;
  352.     }
  353.   fa=cgetg(lx,17);for(i=1;i<lx;i++) fa[i]=zero;
  354.   d=content(x);vv=varn(x);a=gdiv(x,d);ap=deriv(a,vv);t=ggcd(a,ap);v=gdiv(a,t);
  355.   w=gdiv(ap,t);j=0;f=1;nbfac=0;
  356.   while(f)
  357.     {
  358.       j++;w=gsub(w,deriv(v,vv));f=signe(w);
  359.       if(f)
  360.     {
  361.       res=ggcd(v,w);v=gdiv(v,res);w=gdiv(w,res);
  362.     }
  363.       else res=v;
  364.       fa[j]=(lgef(res)>3) ? (long)squff(res,klim) : lgetg(1,18);
  365.       nbfac+=(lg(fa[j])-1);
  366.     }
  367.   av2=avma;y=cgetg(3,19);p1=cgetg(nbfac+1,18);y[1]=(long)p1;
  368.   p2=cgetg(nbfac+1,18);y[2]=(long)p2;
  369.   for(i=1,k=0;i<=j;i++)
  370.     for(i1=1;i1<lg(fa[i]);i1++)
  371.       {
  372.     p1[++k]=lcopy(((GEN)fa[i])[i1]);p2[k]=lstoi(i);
  373.       }
  374.   return gerepile(av,av2,y);
  375. }
  376.  
  377. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  378. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  379. /*~                                                                     ~*/
  380. /*~                            FACTORISATION                            ~*/
  381. /*~                                                                     ~*/
  382. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  383. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  384.  
  385. GEN   factor(x)
  386.   
  387.    GEN   x;
  388.   
  389. {
  390.   long tx=typ(x),l,tetpil,i;
  391.   GEN  y,p1,p2,a0,p,p3,p4,p5;
  392.  
  393.   if(gcmp0(x))
  394.     {
  395.       y=cgetg(3,19);p1=cgetg(2,18);p2=cgetg(2,18);y[1]=(long)p1;
  396.       y[2]=(long)p2;p1[1]=zero;p2[1]=un;return y;
  397.     }
  398.   switch(tx)
  399.   {
  400.   case 1 : y=decomp(x);break;
  401.   case 5 : l=avma;x=gred(x);
  402.   case 4 : if(tx==4) l=avma;p1=decomp(x[1]);
  403.     p2=decomp(x[2]);p4=concat(p1[1],p2[1]);
  404.     p5=concat(p1[2],gneg(p2[2]));p3=indexsort(p4);
  405.     tetpil=avma;y=cgetg(3,19);y[1]=(long)extract(p4,p3);
  406.     y[2]=(long)extract(p5,p3);y=gerepile(l,tetpil,y);break;
  407.   case 10 : p1=content(x);if(!gcmp1(p1)) x=gdiv(x,p1);
  408.     a0=(GEN)(x[2]);
  409.     if(typ(a0)==3)
  410.     {
  411.       p=(GEN)(a0[1]);y=factmod(x,p);
  412.     }
  413.     else
  414.     {
  415.       if(typ(a0)==1) y=factpol(x,0);
  416.       else err(impl,"factor of general polynomial");
  417.     }
  418.     break;
  419.   case 14 : l=avma;x=gred(x);
  420.   case 13 : if(tx==13) l=avma;p1=factor(x[1]);
  421.     p2=factor(x[2]);p3=gneg(p2[2]);tetpil=avma;
  422.     y=cgetg(3,19);y[1]=lconcat(p1[1],p2[1]);
  423.     y[2]=lconcat(p1[2],p3);
  424.     y=gerepile(l,tetpil,y);break;
  425.   case 17:
  426.   case 18:
  427.   case 19: l=lg(x);y=cgetg(l,tx);
  428.     for(i=1;i<l;i++) y[i]=(long)factor(x[i]);break;
  429.   default: err(impl,"general factorization");
  430.   }
  431.   return y;
  432. }
  433.  
  434.  
  435.  
  436. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  437. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  438. /*                                                                 */
  439. /*                         PGCD GENERAL                            */
  440. /*                                                                 */
  441. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  442. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  443.  
  444.  
  445. GEN   ggcd(x,y)
  446.   
  447.    GEN   x,y;
  448.   
  449. {
  450.   long l,l1,tetpil,i;
  451.   long tx=typ(x),ty=typ(y),vx,vy;
  452.   GEN p1,p2,z;
  453.  
  454.   if(tx>ty) {p1=x;x=y;y=p1;l=tx;tx=ty;ty=l;}
  455.   if(ty>=17)
  456.     {l=lg(y);z=cgetg(l,ty);for(i=1;i<l;i++) z[i]=lgcd(x,y[i]);return z;}
  457.   l=avma;
  458.   if(tx<9)
  459.     {
  460.       if(ty<9)
  461.     switch(tx)
  462.       {
  463.       case 1 : switch(ty)
  464.         {
  465.         case 1 : z=mppgcd(x,y);break;
  466.         case 3 : z=cgetg(3,3);
  467.           l=avma;p1=mppgcd(y[1],y[2]);
  468.           if(!gcmp1(p1)) {tetpil=avma;p1=gerepile(l,tetpil,mppgcd(x,p1));}
  469.           z[2]=(long)p1;z[1]=copyifstack(y[1]);
  470.           break;
  471.         case 4 :
  472.         case 5 : z=cgetg(3,4);z[2]=lcopy(y[2]);
  473.           z[1]=lmppgcd(x,y[1]);gredsp(&z);
  474.           break;
  475.         case 7 : z=gpuigs(y[2],min(valp(y),ggval(x,y[2])));break;
  476.         default: z=gun;
  477.         } break;
  478.       case 3 : 
  479.         switch(ty)
  480.           {
  481.           case 3 : 
  482.         z=cgetg(3,3);
  483.         z[1]=gegal(x[1],y[1]) ? copyifstack(x[1]):(long)mppgcd(x[1],y[1]);
  484.         if(gcmp1(z[1])) z[2]=zero;
  485.         else
  486.           {
  487.             l=avma;p1=mppgcd(z[1],x[2]);
  488.             if(!gcmp1(p1))
  489.               {
  490.             tetpil=avma;
  491.             p1=gerepile(l,tetpil,mppgcd(p1,y[2]));
  492.               }
  493.             z[2]=(long)p1;
  494.           } break;
  495.           case 4 :
  496.         p1=mppgcd(x[1],y[2]);
  497.         if(!gcmp1(p1)) err(gcder1);
  498.         tetpil=avma;z=gerepile(l,tetpil,ggcd(y[1],x));
  499.         break;
  500.           case 5 : p1=gred(y);tetpil=avma;
  501.         z=gerepile(l,tetpil,ggcd(p1,z));
  502.         break;
  503.           case 7 : z=gpuigs(y[2],min(valp(y),ggval(x,y[2])));break;
  504.           default: z=gun;
  505.           } break;
  506.       case 4 :
  507.       case 5 : 
  508.         switch(ty)
  509.           {
  510.           case 4:
  511.           case 5: z=cgetg(3,4);
  512.         z[2]=lmulii(x[2],y[2]);l=avma;
  513.         p1=mulii(x[1],y[2]);
  514.         p2=mulii(x[2],y[1]);tetpil=avma;
  515.         z[1]=lpile(l,tetpil,mppgcd(p1,p2));
  516.         gredsp(&z);break;
  517.           case 6:
  518.           case 8: z=gun;break;
  519.           case 7 : z=gpuigs(y[2],min(valp(y),ggval(x,y[2])));break;
  520.           }break;
  521.       case 7: if((ty!=7)||(!gegal(x[2],y[2]))) z=gun;
  522.       else z=gpuigs(y[2],min(valp(y),valp(x)));break;
  523.       default: z=gcmp0(x)?gcopy(y):gun;
  524.       }
  525.       else z=gcmp0(x)?gcopy(y):gun;
  526.     }
  527.   else /* ici tx et ty>=9 */
  528.     {
  529.       vx=gvar9(x);vy=gvar9(y);
  530.       if(vy<vx) return (gcmp0(x))?gcopy(y):polun[vy];
  531.       if(vx<vy) return (gcmp0(y))?gcopy(x):polun[vx];
  532.       switch(tx)
  533.     {
  534.     case 9 : switch(ty)
  535.       {
  536.       case 9 : z=cgetg(3,9);
  537.         z[1]=gegal(x[1],y[1]) ? copyifstack(x[1]):(long)ggcd(x[1],y[1]);
  538.         if(lgef(z[1])<=3) z[2]=zero;
  539.         else
  540.           {
  541.         l=avma;p1=ggcd(z[1],x[2]);
  542.         if(lgef(p1)>3) {tetpil=avma;p1=gerepile(l,tetpil,ggcd(p1,y[2]));}
  543.         z[2]=(long)p1;
  544.           } break;
  545.       case 10: z=cgetg(3,9);z[1]=copyifstack(x[1]);
  546.         l=avma;p1=ggcd(x[1],x[2]);
  547.         if(lgef(p1)>3) {tetpil=avma;p1=gerepile(l,tetpil,ggcd(y,p1));}
  548.         z[2]=(long)p1;
  549.         break;
  550.       case 13:
  551.         p1=ggcd(x[1],y[2]);
  552.         if(!gcmp1(p1)) err(gcder1);
  553.         tetpil=avma;z=gerepile(l,tetpil,ggcd(y[1],x));
  554.         break;
  555.       case 14 : p1=gred(y);tetpil=avma;
  556.         z=gerepile(l,tetpil,ggcd(p1,z));
  557.         break;
  558.       default: err(gcder4);
  559.       } break;
  560.       
  561.     case 10: switch(ty)
  562.       {
  563.       case 10: z=polgcd(x,y);break;
  564.       case 11: z=gpuigs(polx[vx],min(valp(y),gval(x,vx)));
  565.         break;
  566.       case 13: z=cgetg(3,13);z[2]=y[2];z[1]=lgcd(x,y[1]);
  567.         tetpil=avma;z=gerepile(l,tetpil,gred(z));
  568.         break;
  569.       case 14: z=cgetg(3,13);z[2]=lcopy(y[2]);z[1]=lgcd(x,y[1]);break;
  570.       default: err(gcder2);
  571.       } break;
  572.     case 11 : switch(ty)
  573.       {
  574.       case 11: z=gpuigs(polx[vx],min(valp(x),valp(y)));
  575.         break;
  576.       case 13:
  577.       case 14: z=gpuigs(polx[vx],min(valp(x),gval(y,vx)));
  578.         break;
  579.       default: err(gcder2);
  580.       } break;
  581.     case 13 :
  582.     case 14 : if(ty>=15) err(gcder3);
  583.       z=cgetg(3,13);
  584.       z[2]=lmul(x[2],y[2]);l1=avma;
  585.       p1=gmul(x[1],y[2]);
  586.       p2=gmul(x[2],y[1]);tetpil=avma;
  587.       z[1]=lpile(l1,tetpil,ggcd(p1,p2));
  588.       if(ty==13) {tetpil=avma;z=gerepile(l,tetpil,gred(z));}
  589.       break;
  590.     default: err(gcder4);
  591.     }
  592.     }
  593.   return z;
  594. }
  595.  
  596. GEN   glcm(x,y)
  597.      GEN x,y;
  598. {
  599.   long av=avma,tetpil,tx=typ(x),ty=typ(y),i,l;
  600.   GEN p1,p2,z;
  601.  
  602.   if(ty>=17)
  603.     {l=lg(y);z=cgetg(l,typ(y));for(i=1;i<l;i++) z[i]=(long)glcm(x,y[i]);return z;}
  604.   if(tx>=17)
  605.     {l=lg(x);z=cgetg(l,typ(x));for(i=1;i<l;i++) z[i]=(long)glcm(x[i],y);return z;}
  606.   if(gcmp0(x)) return gzero;
  607.   p1=ggcd(x,y);p2=gmul(x,y);
  608.   tetpil=avma;return gerepile(av,tetpil,gdiv(p2,p1));
  609. }
  610.  
  611. GEN   polgcdnun(x,y)
  612.   
  613.    GEN   x,y;
  614.   
  615. {
  616.   long  l,tetpil,tetpil2;
  617.   GEN z,p1,p2,p3;
  618.  
  619.   if(gcmp0(y)) z=gcopy(x);
  620.   else
  621.   {
  622.     p1=x;p2=y;l=avma;tetpil=0;
  623.     while(!gcmp0(p2))
  624.     {
  625.       tetpil2=tetpil;tetpil=avma;p3=gres(p1,p2);
  626.       p1=p2;p2=p3;
  627.     }
  628.     if(tetpil2)
  629.     {
  630.       avma=tetpil;
  631.       if(tetpil2!=l) z=gerepile(l,tetpil2,p1);
  632.       else z=p1;
  633.     }
  634.     else  {avma=l;z=gcopy(y);}
  635.   }
  636.   return z;
  637. }
  638.  
  639. int typecorps(x)
  640.      GEN x;
  641. /* renvoie 1 si probablement un corps simple, 0 sinon */
  642.  
  643. {
  644.   switch(typ(x))
  645.     {
  646.     case 2:
  647.     case 3:
  648.     case 7:
  649.     case 11: return 1;
  650.     case 6: return typecorps(x[1])|typecorps(x[2]);
  651.     default: return 0;
  652.     }
  653. }
  654.  
  655. GEN    polgcd(x,y)
  656.   
  657.    GEN x,y;
  658.   
  659. {
  660.   GEN  p1,p2;
  661.   long l,tetpil,v,e,lx,ty;
  662.  
  663.   if((typ(x)!=10)||(typ(y)!=10)) err(polgcder1);
  664.   if(!signe(y)) return gcopy(x);
  665.   if(!signe(x)) return gcopy(y);
  666.   l=avma;
  667.   if((v=varn(x))==varn(y))
  668.   {
  669.     if(ismonome(x))
  670.     {
  671.       lx=lgef(x);e=gval(y,v);if((lx-3)<e) e=lx-3;
  672.       p1=ggcd(x[lgef(x)-1],content(y));p2=gpuigs(polx[v],e);
  673.       tetpil=avma;return gerepile(l,tetpil,gmul(p1,p2));
  674.     }
  675.     if(ismonome(y))
  676.     {
  677.       lx=lgef(y);e=gval(x,v);if((lx-3)<e) e=lx-3;
  678.       p1=ggcd(y[lgef(y)-1],content(x));p2=gpuigs(polx[v],e);
  679.       tetpil=avma;return gerepile(l,tetpil,gmul(p1,p2));
  680.     }
  681.   }
  682.   if(typecorps(x[lgef(x)-1])||typecorps(y[lgef(y)-1]))
  683.     p1=polgcdnun(x,y);
  684.   else p1=srgcd(x,y);
  685.   if(gcmp0(p1)) return p1;
  686.   if(typ(p1)==10) 
  687.     {
  688.       ty=typ(p1[lgef(p1)-1]);
  689.       if((ty==3)||(ty>5)) return p1;
  690.       if(gsigne(p1[lgef(p1)-1])<0) 
  691.     {tetpil=avma;return gerepile(l,tetpil,gneg(p1));}
  692.       else return p1;
  693.     }
  694.   else
  695.     {
  696.       tetpil=avma;return gerepile(l,tetpil,gmul(polun[v],p1));
  697.     }
  698. }
  699.  
  700.       
  701.       
  702. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  703. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  704. /*~                                                                     ~*/
  705. /*~                           BEZOUT GENERAL                            ~*/
  706. /*~                                                                     ~*/
  707. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  708. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  709.  
  710. GEN   gbezout(x,y,u,v)
  711.    GEN    x,y,*u,*v;
  712. {
  713.   long tx=typ(x),ty=typ(y);
  714.  
  715.   if(ty==1)
  716.   {
  717.     if(tx==1) return bezout(x,y,u,v);
  718.     if(tx==10) {*u=gzero;*v=gdivsg(1,y);return gun;}
  719.     else err(bezoutpoler);
  720.   }
  721.   if(ty==10)
  722.   {
  723.     if(tx==10) return bezoutpol(x,y,u,v);
  724.     if(tx<10) {*u=gdivsg(1,x);*v=gzero;return gun;}
  725.     else err(bezoutpoler);
  726.   }
  727.   if(ty<10)
  728.   {
  729.     if(tx==10) {*u=gzero;*v=gdivsg(1,y);return gun;}
  730.     else err(bezoutpoler);
  731.   }
  732.   err(bezoutpoler);
  733. }
  734.  
  735. GEN   vecbezout(x,y)
  736.    GEN   x,y;
  737. {
  738.   GEN z;
  739.  
  740.   z=cgetg(4,17);z[3]=(long)gbezout(x,y,z+1,z+2);
  741.   return z;
  742. }
  743.  
  744. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  745. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  746. /*                                                                 */
  747. /*                    CONTENU ET PARTIE PRIMITIVE                  */
  748. /*                                                                 */
  749. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  750. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  751.  
  752.  
  753. GEN   content(x)
  754.    GEN   x;
  755.   
  756. {
  757.   long  l,lx,tetpil,i,f,tx=typ(x),txi;
  758.   GEN p1,p2;
  759.  
  760.   if(tx<10)
  761.     {
  762.       if((tx==2)||((tx>=6)&&(tx<=9))) return gun;
  763.       else return gcopy(x);
  764.     }
  765.   lx=lg(x);
  766.   switch(tx)
  767.     {
  768.     case 13:
  769.     case 14: l=avma;p1=content(x[1]);p2=content(x[2]);
  770.       tetpil=avma;return gerepile(l,tetpil,gdiv(p1,p2));
  771.     case 19: if(lx==1) return gun;
  772.       tetpil=l=avma;p1=content(x[1]);
  773.       for(i=2;i<lx;i++) {tetpil=avma;p1=ggcd(p1,content(x[i]));}
  774.       return gerepile(l,tetpil,p1);
  775.     case 10: lx=lgef(x);
  776.     case 11: if((!signe(x))&&((tx==11)||(lx==2))) return gzero;
  777.     default: f=1;
  778.       if(tx!=15) {for(i=lontyp[tx];(i<lx)&&f;i++) f=(typ(x[i])==1);i=lx-1;}
  779.       else i=lx-2;
  780.       p1=(GEN)x[i];l=avma;
  781.       if(f)
  782.     {
  783.       while((i>lontyp[tx])&&(!gcmp1(p1)))
  784.         {--i;tetpil=avma;p1=mppgcd(p1,x[i]);}
  785.     }
  786.       else
  787.     {
  788.       txi=typ(p1);
  789.       if((txi==2)||((txi>=6)&&(txi<=9))) return gun;
  790.       i--;for(;i>=lontyp[tx];i--)
  791.         {
  792.           txi=typ(x[i]);
  793.           if((txi==2)||((txi>=6)&&(txi<=9))) {avma=l;return gun;}
  794.           else {tetpil=avma;p1=ggcd(p1,x[i]);}
  795.         }
  796.     }
  797.       if(l==avma) return gcopy(p1);
  798.       else return gerepile(l,tetpil,p1);
  799.     }
  800. }
  801.  
  802.  
  803. GEN   primpart(x)
  804.   
  805.    GEN   x;
  806.   
  807. {
  808.   long  l,tetpil;
  809.   GEN z,p1;
  810.  
  811.   if(!signe(x)) z=gcopy(x);
  812.   else
  813.   {
  814.     l=avma;p1=content(x);
  815.     if(gcmp1(p1)) {avma=l;z=gcopy(x);}
  816.     else
  817.     {
  818.       tetpil=avma;
  819.       z=gerepile(l,tetpil,gdiv(x,p1));
  820.     }
  821.   }
  822.   return z;
  823. }
  824.  
  825. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  826. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  827. /*                                                                 */
  828. /*                         SOUS RESULTANT                          */
  829. /*                                                                 */
  830. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  831. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  832.  
  833.  
  834. GEN   subres(x,y)
  835.   
  836.    GEN   x,y;
  837.   
  838. {
  839.   long  degq,av,tetpil,f,tx=typ(x),ty=typ(y),dx,dy,du,dv,dr,signh;
  840.   GEN z,g,h,r,p1,p2,p3,p4,u,v;
  841.  
  842.   if(gcmp0(x)||gcmp0(y)) return gzero;
  843.   if((tx<10)||(ty<10))
  844.     {
  845.       if(tx==10) return gpuigs(y,lgef(x)-3);
  846.       if(ty==10) return gpuigs(x,lgef(y)-3);
  847.       else return gun;
  848.     }
  849.   if((tx!=10)||(ty!=10)) err(subrer1);
  850.   if(varn(x)!=varn(y))
  851.     return (varn(x)<varn(y))?gpuigs(y,lgef(x)-3):gpuigs(x,lgef(y)-3);
  852.   dx=lgef(x);dy=lgef(y);
  853.   if (dx<dy)
  854.   {p1=x;x=y;y=p1;f=dx;dx=dy;dy=f;}
  855.   av=avma;p4=content(y);
  856.   if(dy==3) {tetpil=avma;return gerepile(av,tetpil,gpuigs(p4,dx-3));}
  857.   p3=content(x);
  858. /*  p3=gun;p4=gun; */
  859.   u=gdiv(x,p3);v=gdiv(y,p4);
  860.   g=gun;h=gun;f=1;signh=1;
  861.   while (f)
  862.   {
  863.     du=lgef(u);dv=lgef(v);degq=du-dv;
  864.     r=psres(u,v);dr=lgef(r);
  865.     if(dr<=2) f=0;
  866.     else
  867.     {
  868.       u=v;
  869.       if(degq==1) p1=gmul(h,g);
  870.       else p1=gmul(gpuigs(h,degq),g);
  871.       v=gdiv(r,p1);
  872.       g=(GEN)u[lgef(u)-1];
  873.       if(degq==1) h=g;
  874.       else if(degq)
  875.       {
  876.         p1=gpuigs(g,degq);p2=gpuigs(h,degq-1);
  877.         h=gdiv(p1,p2);
  878.       }
  879.       if(((du-3)*(dv-3))&1) signh= -signh;
  880.       if(dr==3) f=0;
  881.     }
  882.   }
  883.   if(dr==2) {z=gzero;avma=av;}
  884.   else
  885.   {
  886.     if(dv==4)
  887.     {
  888.       tetpil=avma;p2=gcopy(v[2]);
  889.     }
  890.     else
  891.     {
  892.       if(dv-3)
  893.       {
  894.         p1=gpuigs(v[2],dv-3);tetpil=avma;
  895.         p2=gdiv(p1,gpuigs(h,dv-4));
  896.       }
  897.       else
  898.       {
  899.         tetpil=avma;p2=gcopy(h);
  900.       }
  901.     }
  902.     if(!gcmp1(p3))
  903.     {
  904.       p1=gpuigs(p3,dy-3);tetpil=avma;p2=gmul(p2,p1);
  905.     }
  906.     if(!gcmp1(p4))
  907.     {
  908.       p1=gpuigs(p4,dx-3);tetpil=avma;p2=gmul(p2,p1);
  909.     }
  910.     if(signh<0) {tetpil=avma;p2=gneg(p2);}
  911.     z=gerepile(av,tetpil,p2);
  912.   }
  913.   return z;
  914. }
  915.  
  916.  
  917. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  918. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  919. /*                                                                 */
  920. /*               RESULTANT PAR MATRICE DE SYLVESTER                */
  921. /*                                                                 */
  922. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  923. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  924.  
  925. GEN resultant2(x,y)
  926.      GEN x,y;
  927. {
  928.   long av=avma,tetpil,dx,dy,i,j,tx=typ(x),ty=typ(y),f;
  929.   GEN p1,p2;
  930.  
  931.  
  932.   if(gcmp0(x)||gcmp0(y)) return gzero;
  933.   if((tx<10)||(ty<10))
  934.     {
  935.       if(tx==10) return gpuigs(y,lgef(x)-3);
  936.       if(ty==10) return gpuigs(x,lgef(y)-3);
  937.       else return gun;
  938.     }
  939.   if((tx!=10)||(ty!=10)) err(subrer1);
  940.   if(varn(x)!=varn(y))
  941.     return (varn(x)<varn(y))?gpuigs(y,lgef(x)-3):gpuigs(x,lgef(y)-3);
  942.   dx=lgef(x)-3;dy=lgef(y)-3;
  943.   if (dx<dy)
  944.   {p1=x;x=y;y=p1;f=dx;dx=dy;dy=f;}
  945.   p1=cgetg(dx+dy+1,19);
  946.   for(j=1;j<=dy;j++)
  947.     {
  948.       p2=cgetg(dx+dy+1,18);p1[j]=(long)p2;
  949.       for(i=1;i<j;i++) p2[i]=zero;
  950.       for(i=j;i<=j+dx;i++) p2[i]=x[dx-i+j+2];
  951.       for(i=j+dx+1;i<=dx+dy;i++) p2[i]=zero;
  952.     }
  953.   for(j=1;j<=dx;j++)
  954.     {
  955.       p2=cgetg(dx+dy+1,18);p1[j+dy]=(long)p2;
  956.       for(i=1;i<j;i++) p2[i]=zero;
  957.       for(i=j;i<=j+dy;i++) p2[i]=y[dy-i+j+2];
  958.       for(i=j+dy+1;i<=dx+dy;i++) p2[i]=zero;
  959.     }
  960.   tetpil=avma;return gerepile(av,tetpil,det(p1));
  961. }
  962.  
  963.  
  964. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  965. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  966. /*                                                                 */
  967. /*           P.G.C.D PAR L'ALGORITHME DU SOUS RESULTANT            */
  968. /*                                                                 */
  969. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  970. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  971.  
  972.  
  973. GEN  srgcd(x,y)
  974.   
  975.    GEN   x,y;
  976.   
  977. {
  978.   long  degq,av,av1,tetpil,f,vx,tx=typ(x),ty=typ(y),dx,dy,du,dv,dr,com;
  979.   GEN d,g,h,r,p1,p2,p3,p4,u,v;
  980.  
  981.   if(gcmp0(x)||gcmp0(y)) return gzero;
  982.   if((tx<10)||(ty<10)) return gun;
  983.   if((tx!=10)||(ty!=10)) err(subrer1);
  984.   if((vx=varn(x))!=varn(y)) return gun;
  985.   dx=lgef(x);dy=lgef(y);
  986.   if (dx<dy)
  987.   {p1=x;x=y;y=p1;f=dx;dx=dy;dy=f;}
  988.   av=avma;
  989.   p3=content(x);p4=content(y);d=ggcd(p3,p4);
  990.   tetpil=avma;d=gmul(d,polun[vx]);
  991.   if(dy==3) return gerepile(av,tetpil,d);
  992.   av1=avma;u=gdiv(x,p3);v=gdiv(y,p4);g=gun;h=gun;f=1;com=0;
  993.   while (f==1)
  994.   {
  995.     du=lgef(u);dv=lgef(v);degq=du-dv;com++;
  996.     r=psres(u,v);dr=lgef(r);
  997.     if(dr<=3) {f=(dr==3)?2:0;}
  998.     else
  999.     {
  1000.       if(com==10)
  1001.       {
  1002.         com=0;u=gdiv(v,content(v));v=gdiv(r,content(r));
  1003.         g=gun;h=gun;
  1004.       }
  1005.       else
  1006.       {
  1007.         u=v;
  1008.         if(degq==1) p1=gmul(h,g);
  1009.         else p1=gmul(gpuigs(h,degq),g);
  1010.         v=gdiv(r,p1);
  1011.         g=(GEN)u[lgef(u)-1];
  1012.         if(degq==1) h=g;
  1013.         else if(degq)
  1014.         {
  1015.           p1=gpuigs(g,degq);p2=gpuigs(h,degq-1);
  1016.           h=gdiv(p1,p2);
  1017.         }
  1018.       }
  1019.     }
  1020.   }
  1021.   if(f) {avma=av1;return gerepile(av,tetpil,d);}
  1022.   p1=gdiv(v,content(v));tetpil=avma;
  1023.   return gerepile(av,tetpil,gmul(d,p1));
  1024. }
  1025.  
  1026.  
  1027. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1028. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1029. /*                                                                 */
  1030. /*                PSEUDO-DIVISION DES POLYNOMES                    */
  1031. /*                                                                 */
  1032. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1033. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1034.  
  1035.  
  1036. GEN   psres(x,y)
  1037.   
  1038.    GEN   x,y;
  1039.   
  1040. {
  1041.   long  l,tetpil,degx,degy;
  1042.   GEN z,p1,p2;
  1043.  
  1044.   degx=lgef(x);degy=lgef(y);
  1045.   if (degx<degy) err(poler15);
  1046.   l=avma;
  1047.   p2=gpuigs(y[degy-1],degx-degy+1);
  1048.   p1=gmul(p2,x);
  1049.   tetpil=avma;
  1050.   z=gerepile(l,tetpil,gres(p1,y));
  1051.   return z;
  1052. }
  1053.  
  1054.  
  1055. GEN   discsr(x)
  1056.   
  1057.    GEN   x;
  1058.   
  1059. {
  1060.   long  dx,av=avma,tetpil;
  1061.   GEN z,p1,p2;
  1062.  
  1063.   switch(typ(x))
  1064.   {
  1065.   case 6: z=stoi(-4);break;
  1066.   case 8: z=discsr(x[1]);break;
  1067.   case 10: dx=lgef(x);p1=deriv(x,varn(x));
  1068.     p1=subres(x,p1);tetpil=avma;p1=gdiv(p1,x[dx-1]);
  1069.     if (((dx-3)&3)>1)
  1070.     {tetpil=avma;z=gerepile(av,tetpil,gneg(p1));}
  1071.     else z=gerepile(av,tetpil,p1);break;
  1072.   case 9: z=discsr(x[1]);break;
  1073.   case 15:
  1074.   case 16: p1=mulii(x[2],x[2]);p2=shifti(mulii(x[1],x[3]),2);tetpil=avma;
  1075.     z=gerepile(av,tetpil,subii(p1,p2));break;
  1076.   default: err(discsrer1);
  1077.   }
  1078.   return z;
  1079. }
  1080.  
  1081. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1082. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1083. /*~                                    ~*/
  1084. /*~                  ALGORITHME DE STURM                               ~*/
  1085. /*~                                    ~*/
  1086. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1087. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1088.  
  1089. long sturm(x)
  1090.   
  1091.    GEN   x;
  1092.   
  1093. {
  1094.   long  degq,av,tx=typ(x),dx,du,dv,dr,sr,s,t,r1,lr;
  1095.   GEN g,h,r,p1,p2,u,v;
  1096.  
  1097.   if(gcmp0(x)||(tx!=10)) err(poltyper);
  1098.   dx=lgef(x);
  1099.   if(dx==3) return 0;
  1100.   av=avma;
  1101.   u=gdiv(x,content(x));v=deriv(x,varn(x));gdiv(v,content(v));
  1102.   g=gun;h=gun;s=gsigne(u[lgef(u)-1]);r1=1;
  1103.   t=(lgef(u)&1) ? -s:s;
  1104.   for(;;)
  1105.   {
  1106.     du=lgef(u);dv=lgef(v);degq=du-dv;
  1107.     r=psres(u,v);dr=lgef(r);
  1108.     if(dr<=2) err(sturmer2);
  1109.     if((signe(v[lgef(v)-1])>0)||(degq&1)) r=gneg(r);
  1110.     lr=lgef(r)-1;
  1111.     sr=gsigne(r[lr]);if(sr!=s) {s= -s;r1--;}
  1112.     if(lr&1) sr= -sr;if(sr!=t) {t= -t;r1++;}
  1113.     if(lr==2) {avma=av;return r1;}
  1114.     u=v;
  1115.     if(degq==1) p1=gmul(h,g);else p1=gmul(gpuigs(h,degq),g);
  1116.     v=gdiv(r,p1);g=gabs(u[lgef(u)-1]);
  1117.     if(degq==1) h=g;
  1118.     else if(degq)
  1119.       {
  1120.         p1=gpuigs(g,degq);p2=gpuigs(h,degq-1);
  1121.         h=gdiv(p1,p2);
  1122.       }
  1123.   }
  1124. }
  1125.  
  1126. long sturmpart(x,a,b)
  1127.   
  1128.    GEN   x,a,b;
  1129.   
  1130. {
  1131.   long  degq,av,tx=typ(x),dx,du,dv,dr,sr,s,t,r1;
  1132.   GEN g,h,r,p1,p2,u,v;
  1133.  
  1134.   if(gcmp0(x)||(tx!=10)) err(poltyper);
  1135.   dx=lgef(x);
  1136.   if(dx==3) return 0;
  1137.   av=avma;
  1138.   u=gdiv(x,content(x));v=deriv(x,varn(x));gdiv(v,content(v));
  1139.   g=gun;h=gun;s=gsigne(poleval(u,b));t=gsigne(poleval(u,a));
  1140.   sr=gsigne(poleval(v,b));
  1141.   r1=0;if(sr!=s) {s= -s;r1--;} 
  1142.   sr=gsigne(poleval(v,a));
  1143.   if(sr!=t) {t= -t;r1++;}
  1144.   for(;;)
  1145.   {
  1146.     du=lgef(u);dv=lgef(v);degq=du-dv;
  1147.     r=psres(u,v);dr=lgef(r);
  1148.     if(dr<=2) err(sturmer2);
  1149.     if((gsigne(v[lgef(v)-1])>0)||(degq&1)) r=gneg(r);
  1150.     sr=gsigne(poleval(r,b));if(sr!=s) {s= -s;r1--;}
  1151.     sr=gsigne(poleval(r,a));if(sr!=t) {t= -t;r1++;}
  1152.     if(dr==3) {avma=av;return r1;}
  1153.     u=v;
  1154.     if(degq==1) p1=gmul(h,g);else p1=gmul(gpuigs(h,degq),g);
  1155.     v=gdiv(r,p1);g=gabs(u[lgef(u)-1]);
  1156.     if(degq==1) h=g;
  1157.     else if(degq)
  1158.       {
  1159.         p1=gpuigs(g,degq);p2=gpuigs(h,degq-1);
  1160.         h=gdiv(p1,p2);
  1161.       }
  1162.   }
  1163. }
  1164.  
  1165. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1166. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1167. /*                                                                 */
  1168. /*         POLYNOME QUADRATIQUE ASSOCIE A UN DISCRIMINANT          */
  1169. /*                                                                 */
  1170. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1171. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1172.  
  1173. GEN   quadpoly(x)
  1174.   
  1175.    GEN   x;
  1176.   
  1177. {
  1178.   long    res,l,tetpil;
  1179.   GEN   y,p1;
  1180.  
  1181.   if(typ(x)!=1) err(arither1);
  1182.   y=cgetg(5,10);
  1183.   y[1]=0x01000005;y[4]=un;l=avma;
  1184.   res=x[lgef(x)-1]&3;
  1185.   if((signe(x)<0)&&res) res=4-res;
  1186.   if(res>1) err(quader);
  1187.   p1=shifti(x,-2);
  1188.   tetpil=avma;
  1189.   if(res)
  1190.   {
  1191.     y[2] = lpile(l,tetpil, (signe(x)<0) ? gsub(un,p1) : gneg(p1));
  1192.     y[3] = lneg(un);
  1193.   }
  1194.   else
  1195.   {
  1196.     y[2] = lpile(l,tetpil,gneg(p1));
  1197.     y[3] = zero;
  1198.   }
  1199.   return y;
  1200. }
  1201.  
  1202. GEN quadgen(x)
  1203.   GEN x;
  1204. {
  1205.   GEN y=cgetg(4,8);y[1]=lquadpoly(x);
  1206.   y[2]=zero;y[3]=un;
  1207.   return y;
  1208. }
  1209.  
  1210. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1211. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1212. /*                                                                 */
  1213. /*                    INVERSE MODULO GENERAL                       */
  1214. /*                                                                 */
  1215. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1216. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1217.  
  1218. GEN   ginvmod(x,y)
  1219.   
  1220.    GEN   x,y;
  1221.   
  1222. {
  1223.   long tx=typ(x),ty=typ(y);
  1224.  
  1225.   if(ty==1)
  1226.   {
  1227.     if(tx==1) return mpinvmod(x,y);
  1228.     if(tx==10) return gzero;
  1229.     else err(ginvmoder);
  1230.   }
  1231.   if(ty==10)
  1232.   {
  1233.     if(tx==10) return polinvmod(x,y);
  1234.     if(tx<10) return gdivsg(1,x);
  1235.     else err(ginvmoder);
  1236.   }
  1237.   err(ginvmoder);
  1238. }
  1239.  
  1240. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1241. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1242. /*~                                    ~*/
  1243. /*~                  POLYGONE DE NEWTON                ~*/
  1244. /*~                                    ~*/
  1245. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1246. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1247.  
  1248. GEN newtonpoly(x,p)
  1249.      GEN x,p;
  1250. {
  1251.   GEN y;
  1252.   long n,*vval,i,a,b,c,u1,u2,r1,r2;
  1253.  
  1254.   if(typ(x)!=10) err(newter1);
  1255.   n=lgef(x)-3;if(n<=0) {y=cgetg(1,17);return y;}
  1256.   vval=(long *)newbloc(4*(n+1));
  1257.   for(i=0;i<=n;i++) vval[i]=(gcmp0(x[i+2])) ? LARGERINT : ggval(x[i+2],p);
  1258.   a=0;b=1;y=cgetg(n+1,17);
  1259.   while(b<=n)
  1260.     {
  1261.       u1=vval[a]-vval[b];u2=b-a;
  1262.       for(c=b+1;c<=n;c++)
  1263.     {
  1264.       r1=vval[a]-vval[c];r2=c-a;
  1265.       if(u1*r2<=u2*r1) {u1=r1;u2=r2;b=c;}
  1266.     }
  1267.       for(i=a+1;i<=b;i++) y[i]=ldiv(stoi(u1),stoi(u2));
  1268.       a=b;b=a+1;
  1269.     }
  1270.   killbloc(vval);return y;
  1271. }
  1272.  
  1273.  
  1274. GEN oldpolfnf(a,t)
  1275.      GEN a,t;
  1276.  
  1277. /* Factorisation du polynome a sur le corps de nombres defini par le
  1278. polynome t */
  1279. {
  1280.   GEN y,p1,p2,u,g,fa,n,r,unt,f,b,pro;
  1281.   long av=avma,tetpil,lx,v,i,e,k;
  1282.   
  1283.   if((typ(a)!=10)||(typ(t)!=10)) err(polfnfer1);
  1284.   if(varn(t)) err(polfnfer2);
  1285.   if(gcmp0(a)) return gcopy(a);
  1286.   v=varn(a);unt=gmodulcp(un,t);
  1287.   u=gdiv(a,ggcd(a,deriv(a,v)));u=gmul(unt,u);
  1288.   setvarn(u,255);g=lift(u);setvarn(u,0);k= -2;
  1289.   do {k++;n=subres(t,gsubst(g,255,gsub(polx[255],gmulsg(k,polx[0]))));}
  1290.   while(!issquarefree(n));
  1291.   fa=(GEN)(factor(n)[1]);lx=lg(fa);y=cgetg(3,19);p1=cgetg(lx,18);y[1]=(long)p1;
  1292.   p2=cgetg(lx,18);y[2]=(long)p2;setvarn(a,0);
  1293.   for(i=1;i<lx;i++)
  1294.     {
  1295.       setvarn(fa[i],0);
  1296.       f=gsubst(fa[i],0,gadd(polx[0],gmulsg(k,gmodulcp(polx[0],t))));
  1297.       pro=ggcd(u,gmul(unt,f));p1[i]=ldiv(pro,pro[lgef(pro)-1]);
  1298.       e=0;b=poldivres(a,p1[i],&r);
  1299.       while(gcmp0(r)) {a=b;e++;b=poldivres(a,p1[i],&r);}
  1300.       p2[i]=lstoi(e);
  1301.     }
  1302.   tetpil=avma;setvarn(a,v);return gerepile(av,tetpil,gcopy(y));
  1303. }
  1304.  
  1305. GEN polfnf(a,t)
  1306.      GEN a,t;
  1307.  
  1308. /* Factorisation du polynome a sur le corps de nombres defini par le
  1309. polynome t */
  1310. {
  1311.   GEN y,p1,p2,u,g,fa,n,r,unt,f,b,pro,ain;
  1312.   long av=avma,tetpil,lx,v,i,e,k;
  1313.   
  1314.   if((typ(a)!=10)||(typ(t)!=10)) err(polfnfer1);
  1315.   if(varn(t)) err(polfnfer2);
  1316.   if(gcmp0(a)) return gcopy(a);
  1317.   v=varn(a);ain=a;setvarn(a,0);unt=gmodulcp(un,t);
  1318.   u=gdiv(a,ggcd(a,deriv(a,0)));u=gmul(unt,u);
  1319.   setvarn(u,255);g=lift(u);setvarn(u,0);k= -2;
  1320.   do {k++;n=subres(t,gsubst(g,255,gsub(polx[255],gmulsg(k,polx[0]))));}
  1321.   while(!issquarefree(n));
  1322.   fa=(GEN)(factor(n)[1]);lx=lg(fa);y=cgetg(3,19);p1=cgetg(lx,18);y[1]=(long)p1;
  1323.   p2=cgetg(lx,18);y[2]=(long)p2;
  1324.   for(i=1;i<lx;i++)
  1325.     {
  1326.       setvarn(fa[i],0);
  1327.       f=gsubst(fa[i],0,gadd(polx[0],gmulsg(k,gmodulcp(polx[0],t))));
  1328.       pro=ggcd(u,gmul(unt,f));
  1329.       p1[i]=(typ(pro)==10)?ldiv(pro,pro[lgef(pro)-1]):(long)pro;
  1330.       e=0;b=poldivres(a,p1[i],&r);
  1331.       while(gcmp0(r)) {a=b;e++;b=poldivres(a,p1[i],&r);}
  1332.       p2[i]=lstoi(e);if(v) p1[i]=lsubst(p1[i],0,polx[v]);
  1333.     }
  1334.   tetpil=avma;setvarn(ain,v);return gerepile(av,tetpil,gcopy(y));
  1335. }
  1336.  
  1337. GEN nfiso(a,b)
  1338.      GEN a,b;
  1339. {
  1340.   long av=avma,tetpil,n,m,i,c,va,vb,lx;
  1341.   GEN p1,y,ain,bin,la,lb,da,db,fa,ex;
  1342.  
  1343.   if((typ(a)!=10)||(typ(b)!=10)) err(nfisoer1);
  1344.   n=lgef(a)-3;m=lgef(b)-3;if((n<=0)||(m<=0)) err(nfisoer2);
  1345.   if(n!=m) return gzero;
  1346.   ain=a;bin=b;va=varn(a);vb=varn(b);setvarn(a,0);setvarn(b,0);
  1347.   p1=content(a);if(!gcmp1(p1)) a=gdiv(a,content(a));
  1348.   p1=content(b);if(!gcmp1(p1)) b=gdiv(b,content(b));
  1349.   la=(GEN)a[n+2];
  1350.   if(!gcmp1(la)) a=gmul(gpuigs(la,n-1),gsubst(a,0,gdiv(polx[0],la)));
  1351.   lb=(GEN)b[n+2];
  1352.   if(!gcmp1(lb)) b=gmul(gpuigs(lb,n-1),gsubst(b,0,gdiv(polx[0],lb)));
  1353.   da=discsr(a);db=discsr(b);p1=gdiv(da,db);
  1354.   if(typ(p1)==4) p1=gmul(p1[1],p1[2]);
  1355.   if(typ(p1)!=1) err(nfisoer3);
  1356.   if(!carreparfait(p1)) return gzero;
  1357.   fa=polfnf(a,b);ex=(GEN)fa[2];p1=(GEN)fa[1];lx=lg(p1);c=0;
  1358.   for(i=1;i<lx;i++) 
  1359.     {
  1360.       if(!gcmp1(ex[i])) err(nfisoer4);
  1361.       if(lgef(p1[i])==4) c++;
  1362.     }
  1363.   if(!c) return gzero;
  1364.   y=cgetg(c+1,17);c=0;
  1365.   for(i=1;i<lx;i++) if(lgef(p1[i])==4) y[++c]=lneg(lift(((GEN)p1[i])[2]));
  1366.   setvarn(ain,va);setvarn(bin,vb);
  1367.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1368. }
  1369.   
  1370. GEN nfincl(a,b)
  1371.      GEN a,b;
  1372. {
  1373.   long av=avma,tetpil,n,m,i,c,va,vb,lx,q;
  1374.   GEN p1,y,ain,bin,la,lb,da,db,fa,ex;
  1375.  
  1376.   if((typ(a)!=10)||(typ(b)!=10)) err(nfisoer1);
  1377.   m=lgef(a)-3;n=lgef(b)-3;if((n<=0)||(m<=0)) err(nfisoer2);
  1378.   if(n%m) return gzero;
  1379.   ain=a;bin=b;va=varn(a);vb=varn(b);setvarn(a,0);setvarn(b,0);
  1380.   p1=content(a);if(!gcmp1(p1)) a=gdiv(a,content(a));
  1381.   p1=content(b);if(!gcmp1(p1)) b=gdiv(b,content(b));
  1382.   la=(GEN)a[m+2];
  1383.   if(!gcmp1(la)) a=gmul(gpuigs(la,m-1),gsubst(a,0,gdiv(polx[0],la)));
  1384.   lb=(GEN)b[n+2];
  1385.   if(!gcmp1(lb)) b=gmul(gpuigs(lb,n-1),gsubst(b,0,gdiv(polx[0],lb)));
  1386.   q=n/m;
  1387.   da=discsr(a);db=discsr(b);
  1388.   if((typ(da)!=1)||(typ(db)!=1)) err(nfisoer3);
  1389.   fa=factor(da);ex=(GEN)fa[2];p1=(GEN)fa[1];lx=lg(p1);
  1390.   for(i=1;i<lx;i++) 
  1391.     if((itos(ex[i])&1)&&(!divise(db,gpuigs(p1[i],q)))) return gzero;
  1392.   fa=polfnf(a,b);ex=(GEN)fa[2];p1=(GEN)fa[1];lx=lg(p1);c=0;
  1393.   for(i=1;i<lx;i++) 
  1394.     {
  1395.       if(!gcmp1(ex[i])) err(nfisoer4);
  1396.       if(lgef(p1[i])==4) c++;
  1397.     }
  1398.   if(!c) return gzero;
  1399.   y=cgetg(c+1,17);c=0;
  1400.   for(i=1;i<lx;i++) if(lgef(p1[i])==4) y[++c]=lneg(lift(((GEN)p1[i])[2]));
  1401.   setvarn(ain,va);setvarn(bin,vb);
  1402.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1403. }
  1404.   
  1405. long rootsof1(x)
  1406.      GEN x;
  1407. {
  1408.   long av=avma,n,c,i,j,k,l,lx,fl,fl2,r1,p=0,pj,w,*lis;
  1409.   byteptr pt=diffptr;
  1410.   GEN dk,fa;
  1411.  
  1412.   if(typ(x)!=10) err(rootofer1);
  1413.   n=lgef(x)-3;if(n<=0) err(rootofer2);
  1414.   r1=sturm(x);if(r1) {avma=av;return 2;}
  1415.   dk=discf(x);w=1;lis=(long*)malloc(100);c=0;
  1416.   p+=*pt++;
  1417.   do
  1418.     {
  1419.       if(!(n%(p-1))) lis[++c]=p;
  1420.       p+=*pt++;
  1421.     }
  1422.   while(p-1<=n);
  1423.   for(i=1;i<=c;i++)
  1424.     {
  1425.       p=lis[i];k=(ggval(dk,stoi(p))+(n/(p-1)))/n;fl=1;
  1426.       if(p==2) {pj=4;j=2;} else {pj=p;j=1;}
  1427.       do
  1428.     {
  1429.       fa=(GEN)polfnf(cyclo(pj),x)[1];lx=lg(fa);fl2=1;
  1430.       for(l=1;(l<lx)&&fl2;l++) fl2=(lgef(fa[l])!=4);
  1431.       if(fl2) {w*=(pj/p);fl=0;}
  1432.       j++;pj*=p;
  1433.     }
  1434.       while((j<=k)&&fl);
  1435.       if(fl) w*=(pj/p);
  1436.     }
  1437.   avma=av;free(lis);return w;
  1438. }
  1439.